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We apply a spin-coherent states formalism to study the central-spin model with monochromatic 
bath and symmetric coupling (the Mermin model); in particular, we derive analytic expressions for 
the integrated density of states in the thermodynamic limit when the number of bath spins is taken 
to infinity. From the thermodynamic limit spectrum we show the phase diagram for the system can 
be divided into four regions, partitioned on the one hand into a symmetric (non-degenerate) phase 
or a broken symmetry (degenerate) phase, and on the other hand by the case of overlapping or 
non-overlapping energy surfaces. The nature and position of singularities appearing in the energy 
surfaces change as one moves from region to region. Our spin-coherent states formalism naturally 
leads us to the Majorana representation, which is useful to transform the Schrodinger equation into 
a Ricatti-like form that can be solved in the thermodynamic limit to obtain closed-form expressions 
for the integrated density of states. The energy surface singularities correspond with critical points 
in the density of states. We then use our results to compute expectation values for the system that 
help to characterize the nature of the quantum phase transition between the symmetric and broken 
phases. 



I. INTRODUCTION 



The study of a two-level system interacting with some form of environment or background has been the focus of 
much research [TH3]. Such studies have focused on important questions such as the onset of decoherence in nature as 
well as the transition from quantum to classical physics. The former question takes on a concrete practical importance 
in the context of the emerging field of quantum information processing, in which case the two-level system could be 
used to simulate the time-evolution of a qubit 4 . In this particular situation, the goal would be to maintain the 
coherence of the two-level system on some sufficiently long time scale in order to preserve the information embedded 
in the qubit until the next step in the computation process occurs. 

While these studies have been considered in the context of a wide range of natural phenomena, interestingly it seems 
these models fall into two distinct classes based on the characterization of the environment [1]. On the one hand, our 
physical picture could be seen as an environment composed of localized impurities, each with a small number of low 
energy states that can be modeled as a spin 1/2 system or a set of such spins. In the case of the two-level central 
system coupled to such a 'spin bath,' this is called the central-spin model. A special case of this model will be the 
focus of the present study. On the other hand we could also consider a two-level system coupled to a large number of 
background oscillators, which is the spin-boson model. 

Here we focus on the so-called Mermin model, which is the special case of the central-spin model in which the 
background spins share a common frequency (the 'monochromatic bath') and each interact symmetrically with the 
central spin. This model was proposed by N. D. Mermin [S] as an illustration of the transition from a delocalized 
quantum state to a localized state (with broken spatial symmetry) in the classical limit. Since then, the dynamical 
properties of the system across the phase transition have been studied [51 [7] , and the model has also been studied 
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under the names of the finite Jaynes-Cummings Model (or su(2) Jaynes-Cummings model) "5] and the spin star system 
[9]. In the latter case, the authors have applied various approximation techniques to the density operator evolution 
in order to demonstrate decoherence and the presence of non-Markovian effects. 

In the present study, we will apply a spin coherent states pi3HT^ based approach to the Mermin model in order to 
obtain the analytic form of the integrated density of states in each phase region of the model. By introducing the 
spin coherent states, we will find that the wave function for the system is naturally transformed into the Majorana 
representation 13J. This approach follows that which was recently applied to analyze the spectrum of the well-known 
Lipkin-Meshkov-Glick (LGM) model [131 [TS]. However, it is worth noting that this type of analysis has proved fruitful 
in other studies as well. For example, a similar approach (also based on the spin coherent states) has previously been 
applied to the Mermin (or finite JCM) model |8], although with the aim of characterizing the system in a somewhat 
different manner, as we will briefly elaborate on in our concluding remarks. As another example, a generalized SU(A/) 
coherent states approach has been recently applied in investigations of the M-site Bose Hubbard model |16j . 

The organization of our paper is as follows. In Sec. [lljwe write the Mermin model Hamiltonian and introduce the 
system representation in terms of the spin coherent states that transforms the wave function into the form of the 
Majorana polynomial. The roots of this polynomial can then be used to characterize the system. In Sec. |III[ we 
apply a mean field theory analysis to find the classical energy sheets in the thermodynamic limit. There are two such 
sheets, one each associated with the two available spin states of the central system. We will further study the extrema 
of these energy sheets and find that their behavior serve to characterize the phase diagram. It is well known that, due 
to spontaneous symmetry breaking, the phase diagram for the Mermin model can be divided into a symmetric phase 
(with non-degenerate ground state) and a (degenerate) broken symmetry phase. We will see that these two regions 
can be further divided according to whether the two energy sheets are overlapping or non-overlapping. 

We will then turn to the analysis of the time- independent Schrodinger equation in Sec. |IV| Using the coherent states 
formalism we show that the Schrodinger equation can be mapped into a first-order nonlinear differential equation 
in a Ricatti-like form. By expanding these equations in powers of the inverse number of background spins 1/iV, 
we formulate a method by which we can determine an analytic form for the integrated density of states in the 
thermodynamic limit (infinite N). The results of these calculations are presented in Sec. [V{ w ith some details of how 
the density of states expressions may be obtained provided in the Appendix. Then in Sec. |VI| we will make use of our 
results to compute the expectation values for several key quantities in the system. In part, this will serve to help in 
the characterization of the quantum phase transition. Finally in Sec. |VII| we summarize our results and make some 
comments on future work. 



II. HAMILTONIAN AND MAJORANA REPRESENTATION 



The Mermin model describes the interaction of a two-level system, represented by the central spin or 'small spin' cr, 
with a collection of N background spins cr*. Considering only the fully symmetric subspace, these background spins 
can be summed into one large 'environment spin' S = X^iLi ""V^ with maximal spin S — N/2. Here the u's are the 
standard Pauli spin matrices, although we will introduce a different representation for the central spin below. The 
interaction term will then be given by an anti-ferromagnetic coupling between the x components of the central spin 
and the environment spin, according to 

/ \ 

^ = T"^+^2^+^^"^25 

with 7x > 0. Also notice that the second and third terms involving the collective spin have been re-scaled by the 
factor 1 /2S in order to obtain a non-trivial phase diagram. 

Let us now introduce the spin coherent states basis that will lead us to the Majorana representation. First, 
we introduce the most familiar representation for this system in terms of the tensor product \S,M) \s,m) of the 
environment spin basis \S, M) and the small spin basis |s, m) by writing a general state as j^*) = ^ CM,m{\S, M) (E) 
|s, m)). Here M = —S, . . . , S and m = — s, . . . , s where, for the system we have in mind, S* 3> s. Then the coherent 
states representation |a) in terms of the collective spin S is given by 

(a|*) =^^CM,™(a|5,M)|s,m) (2) 

m M 

where a is a complex parameter; we now define 



V'™(a) = ^CM.™(a|5,M) =Cn(«-4'"^ (3) 



a — a 

M k=l 
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as the Majorana polynomial of order d < N with maximum order N = 25*. The coherent state \a) is given by 



\a) 



— 5'), which can be shown to be equivalent to 



(25)! 



,3+M 



(S + My.{S - M)l 



S,M). 



The inner product of two spin coherent states is given by 

{a'\a) = + 

In the present representation, the environmental spin operators appearing above are given by 

Sz = —S + ada 

s- = a„ 

S+ = 2Sa~a^da 



(4) 



(5) 



(6) 



along with the ordinary relations S'i = 5*3; ± iSy. 

Now we will choose to specify that s = 1/2 (Mermin's model). Let us write our vector- valued wave function in the 
explicit Pauli representation as 



(n\ - ( '^i/2(a) \^( V't(«) 



(7) 



with the functions il'^,^{a) given through Eq. ([s]). Before proceeding with our analysis, we find that a change of basis 
for the small spin following the unitary transformation 



U = 



1 1 

1 -1 



(8) 



will help simplify the analysis of the time-independent Schrodinger equation in Sec. |IV[ Applying this transformation 
to Eq. ([?]) we find that the wave vector in the 'up, down' representation is transformed into the 'symmetric, anti- 
symmetric' representation according to 



V't(a) 


1 


'0t(") 


+ '0i(a) 




(t>s{a) 




"71 


_ '0t(") 









(9) 



Note that in this basis, the ax = U'^OxU matrix (which includes the interaction term in Eq. ([l]) above) is diagonalized. 
As a result, the Hamiltonian ([!]) takes the form 



H 



23 



7: 



a: 2S 



(10) 



in which operator terms appear only on diagonal entries, which is an indication that the Schrodinger equation may 
be written in a more compact form in this basis. 



III. CLASSICAL ENERGY SURFACES AND PHASE DIAGRAM 



In this section we will introduce a mean field theory (MET) approximation to find the two classical energy surfaces, 
one each associated with the two degrees of freedom for the central spin. From the energy surfaces we can then 
determine the critical points in the spectrum, including the ground state. These are singular points such as saddle 
points that may occur inside the energy continuum or the max/min that form the edges of the continuum for a given 
sheet. The ground state is given by the minimum on the lower energy surface. The phase diagram will then be 
revealed by studying the critical points as a function of the Hamiltonian parameters. 



A. Mean Field Theory and Classical Energy Surfaces 



In the context of the spin coherent states representation we introduce a mean field theory approximation, such that 
correlations between the bath and the central spin factorize, as ^ ~ ^mft{<^) = vi^ + ctct)'^^ with ij = {r]s,r]a) 
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a constant complex vector. With this factorized form we are in essence assuming that in the MFT each individual 
background spin will have its own distinguishable interaction with the central spin, and that these interactions will 
be symmetric. We will further take the thermodynamic limit iV = 25* — > oo to obtain our expressions for the classical 
energy surfaces below. 

Let us make a brief comment on our approximation before proceeding. Note that the thermodynamic limit in this 
context is equivalent to the classical limit in the sense that as 5* — )• oo the large environment spin will behave as a 
classical spin, since we have chosen to work in the symmetric subspace where the individual bath spins can be summed 
into a single giant spin. 

Applying our mean field ansatz, we can obtain the classical energy surfaces by further taking the limit = 25" — > oo 

as 

, + aa) ± Juj'^il + aaf + -^l[a. + 5)2 

2(1 + aa) 

This quantity gives a mean field description of the energy spectrum of our model parameterized as a function of the 
complex pair a, a. The singular points in the spectrum, such as edge singularities and saddle-points, can be exactly 
determined as the critical points of t±{a^ a), a calculation which is carried out in the following section. We plot the 
classical energy surfaces in terms of a; = Re a, y = Im a for four choices of the system parameters in Fig. [2j 

Note that alternatively, we could apply a classical spin description that is equally valid in the thermodynamic limit. 
This classical spin description is obtained via the transformation 

N N N 

Sx ^ — sin 9 cos cjy^Sy — — sin 6 sin cjy^Sz — — cos 9. (12) 

Indeed, the quantum phase transition between the symmetric and broken phases has previously been studied in 
terms of this equivalent description [5HZ]- A transformation between the two representations can be accomplished 
through a — tan |e"^. There are a handful of calculations that are in fact more simple in one representation or the 
other, although we will focus on the a, a description below. To obtain a better physical intuition for these complex 
parameters, we note the following relationship to the classical spin components 

(Re(o))' + (Im(o))^ = ''^J,f^^^f^^ . (15) 



obtained from the above transformation and Eq. (12). Hence, roughly speaking, the real part of a can be associated 



with the X component of spin and the imaginary part with the y component. 

B. Domain of critical points in the o, a plane 

In order to characterize the system and obtain the phase diagram, we determine the critical points of the spectrum 
from da€± — da^± — 0. The solution to these conditions are obtained as roots of a quadratic equation in a^, along 
with the condition a = a (that is, a real). The four explicit solutions of the quadratic are given by ±a_|_,±Q;_, in 
which 



7.^ + »2^2 + 272^2 ± 27.» V(»2 + ^2) ( 2 + ^2) 

= V "ff^^ • ^ ^ 

There are two other sets of critical points given by a = a and a = a — > 00. 

We notice immediately that the denominator 7^ — fi^cj^ = (^y^ _ ftuj2){-y'^ + HlOz) of the critical points given in Eq. 



( 16 1 will vanish whenever 

7' - l^^ll^.l (17) 

is satisfied. Since fl and uj^ may take either positive or negative values, this condition defines four symmetric curves 
in parameter space as shown in Fig. [l] Since the four quadrants given in this diagram are identical, we will focus 
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FIG. 1: The four quadrants in parameter space for the Mermin model. Here we have chosen the value 7a; = 1. We show 

2 

the line \u)z\ ~ ^ that denotes the quantum phase transition in each quadrant. In each quadrant, we have also labeled the 
symmetric (S) and broken phase (BP) portion of the phase diagram. Each quadrant is equivalent according to the symmetry 
of the system, thus we will focus exclusively on Quadrant I as we proceed. Note that we then further divide each quadrant into 
four regions in Fig. [5] below. 



throughout the paper on quadrant I where the above condition can simply be written 7^ — flcuz ■ We wiU demonstrate 
below that this divergence indicates the presence of a quantum phase transition. 

On one side of this transition (w^ < 7^/J7), these four critical points wiU all be real, with ±a_ forming two 
degenerate (global) minima on the lower sheet, while the two points ±a+ will form degenerate (global) maxima in the 
upper sheet. On the other side of the transition {utz > lx/^)t the points a± become complex conjugates (violating 
a = a) and there will be a single minimum (given by a = a = 0) at the origin in the lower sheet, which will form the 
non-degenerate ground state in the symmetric phase. 

Finally, notice from Eq. ( 13 ) that the purely real points in Eq. ( 16 1 are associated with a non-zero x component to 
the spin. Hence, the quantum phase transition is apparently induced by the coupling term in the Hamiltonian (which 
involves the x component of the spins). Writing the condition for the broken phase as 7^ > fluiz makes clear that 
the symmetry breaking can be viewed as occurring when the energy scale of the coupling dominates over that of the 
spins. 



C. Extremal energy values, ground state degeneracy and quantum phase transition 



Having obtained the critical values of a, a, we will now find the explicit values of the extremal energies including 
the ground state. 

We note that the critical points at the origin a — a — are critical points in both sheets throughout parameter 
space. From Eq. (11) we immediately find 

£±(a = 0,a = 0) = ^^^^4 (18) 

The negative choice on the lower sheet acts as the (unique) ground state throughout the symmetric phase. This 
point becomes a saddle point in the broken phase as we show below. Meanwhile, for the pair of critical points at 
infinity a, a ^ co we find 

e-|-(a -» 00, a 00) = ^^'^^ = e±- (19) 



FIG. 2: (a) Region 1 (symmetric, non-overlapping), (b) Region 2 (symmetric, overlapping), (c) Region 3 (broken phase, non- 
overlapping), and (d) Region 4 (broken phase, overlapping). Here we plot the two energy surfaces as a function of j: = Re(a) 
and y = Im(a) for each of the four phase regions in quadrant I. We have used the values (a) Region 1: f2 = 2.2, uj^ = 2.5, (b) 
Region 2: 57 = 2.0, uJz = 0.7, (c) Region Q, = 0.25, Uz — 0.6, and (d) Region 4: = 0.8, ujz = 0.1 along with the typical 
value 7a; = 1 in all cases. Critical points are shown as red points and some level energy surfaces are drawn in black. 



To verify the presence of the quantum phase transition, we construct the Hessian matrix that describes the quah- 
tative behavior of a given critical point through the second derivatives of the classical energy surfaces. We construct 
the Hessian &±{x, y) with x = Rc a,y = Im a according to 



Q±ix,y) 



dx.xe±{x,y) dx,ye±{x,y) 
dy,xf-±{x,y) dy^ye±{x,y) 



(20) 



We can now distinguish between maxima, minima and saddle points, which will have negative definite, positive 
definite, and indefinite Hessian, respectively. 

For the energy value on the negative sheet we calculate the Hessian as 



e_(o,o) = 



2(1] -7^M) 

2n 



(21) 



The determinant is given as det8_(0, 0) = 417(^217 — j'^/uJz)/uJz, which is positive for ojz > 1x1^- Both diagonal 
entries in Eq. (|2ll) are positive for this case, meaning the Hessian is positive definite and the point is indeed a 
minimum. Visual inspection of the energy surface plots in Fig. [2ja, b) reveals this energy to be the unique global 
minimum in the symmetric phase. 

For ujz < the determinant det8_(0, 0) becomes negative, giving an indefinite Hessian so that is now a 

saddle point as can be seen in Figs, ^c, d). This point goes from a stable critical point in the symmetric phase to 
an unstable critical point in the broken phase. We also see in these figures that in this phase a new, twice degenerate 
ground state forms at the points zLa- as we discussed above. These are the stable critical points in the broken phase. 
In parallel with this transition, on the upper sheet the point ej^ goes from a global maximum to a saddle point at 
infinity in the broken phase. 

We can obtain the explicit energy values associated with the broken phase global max/min points by plugging Eq. 



(161 into Eq. (Ill to find 



HP 



^±^Vi^'z+l'xW+l'x)- 
^Ix 



(22) 



The value on the lower sheet e_ gives the energy of the twice degenerate ground state. 



Plugging the critical values for a, a into Eqs. we find the MFT values for the background spin m the 

ground state in either phase. For the symmetric case with a = a = we find {Sx) = (Sy) = while (Sz) ~ N/2 ~ S. 
Meanwhile we can also find the 2 component eigenspinor associated with the symmetric phase ground state eigenvalue 
e° = —{n + uJz)/2 in order to calculate {(j^) = (ay) = and (cr^) = 1. Hence we see that in the symmetric phase 
ground state, the central-spin and the background spin are both aligned along the positive z-axis, as represented in 
the qualitative diagram in Fig. [sj^a). 

Meanwhile for the degenerate ground states in the broken phase we plug Eq. ( 16 1 into Eq. ( 13 1 to find 




(23) 
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FIG. 3: (Color online) Qualitative descriptive diagrams for the ground state in (a) the symmetric phase, and (b) the broken 
phase. The environment spin (large blue arrow) and central spin (black arrow) reside in the x-z plane (dark shading) in either 
case. In the symmetric phase the arrows are aligned along the z axis with the coupling effectively off. In the broken phase, 
two degenerate ground states form with the components of the spin along the a;-axis anti-aligned (the anti-alignment is only 
approximate for the case ojz / f^). 



for the x-component of the background spin as reported in [SHI] ■ We see that the two degenerate ground states are 
associated with a left-right symmetry in terms of this component of the spin. As this component was zero in the 
symmetric phase we can clearly interpret this as a second order quantum phase transition in which plays the role 
of the order parameter. Indeed, we immediately read off the critical exponent for the order parameter from Eq. (23) 
as /3 = 1/2, in agreement with the well-known MFT result [T7]. Here f3 describes the power-law scaling behavior of 
the order parameter in the vicinity of the QPT. 
In a similar manner wc may obtain 



(Sy) 



(24) 



for the y and z components for the degenerate ground state. We note from comparison of Eq. (22) and Eq. (24) that 



the ground state energy in the broken phase is proportional to the mean value of the z component of the background 
spin, according to e^^ = -(w^ + jl){Sz)/2nS. 

Finally we may obtain the expectation values for the small spin as 



= T\ 



= 0, 




(25) 



Comparing the expression for {gx) with {Sx)/S from Eq. (23) we see that the central-spin and the background spin 
in the degenerate ground state are exactly anti- aligned a long the x-axis for the case Wz = with equal energies. For 
this special case we have {Sx)/S = —{(Jx) = ±1/7^ — Vfi /^x (anti-alignment) and {Sz)/S = {(j^) — Q./^x (alignment). 
For the case ^ the exact anti-alignment along x and alignment along z is broken, but it approximately holds for 
the case ~ 17, which is our primary consideration here (we consider that all numerics presented in this paper fall 
into this parameter range). 

The descriptive diagrams in Figs.j3[a, b) may aid us in a qualitative understanding of the behavior of the ground 
state through this transition. Fig. pfa) qualitatively represents the ground state in the symmetric phase with the 
central spin (black arrow) and environment spin (blue arrow) aligned in the positive z axis. Neither spin has an x-axis 
component, which is the coupling axis in our Hamiltonian Eq. ([l]). 

As we increase the relative strength of the coupling ^x and cross into the broken phase we see in Fig. |3]^b) that 
the two spins each acquire a component that lies along the x direction (this diagram is drawn as a representation of 
the exact case uj^ — For example in the left sphere the environment spin has acquired a component along the 
positive X-axis (consistent with the -I- sign choice in Eq. ( |23[ )) while the central spin has acquired a component along 
the negative x-axis (consistent with the — sign choice in Eq. (25)). The spins are anti-aligned in the x direction as 
a result of the anti-ferromagnetic coupling (7^ > 0) in the Hamiltonian (jl) being 'switched on' in the broken phase 
ground state. Of course, there is also the equal-energy configuration with the central spin having a component along 
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r^= <a>^ for |*o(^^=A,n=A,7=l)) 




2.0 



FIG. 4: Behavior of in the thermodynamic limit as a function oi \ = Q = uj^ in the case ~ 1 with numerics provided for 
the S = 10, 20 and 50 cases. 



the positive a:-axis and the environment spin pointing in the negative a;-axis as depicted by the sphere on the right in 
Fig. [Sj^b), which is the source of the broken phase degeneracy. 

Finally, we are in a position to calculate the spin polarization r^, the summed square of the components of the 
spin matrices expectation values (note that for the environment spin, we normalize this quantity by dividing by S). 
This quantity is also known as the one-tangle. Let us describe its properties using the central spin as an example: 
in the case that the central spin is in a pure state, we shall have = 1. In this case, is equivalent to the usual 
Bloch sphere radius. However, to extend the definition of this quantity to the mixed state case, one must allow for 

< 1, in which situation we can uniquely describe the system properties in terms of the Bloch Ball [20]. Then < 1 
serves as a precise measure of the entanglement present in the system; the smaller the value of the greater is the 
entanglement in the system, with r — yielding the completely mixed state [18 . 

For both the central spin and the environment spin, sums to one in the symmetric phase so there is no entangle- 
ment here. However, for the degenerate phase ground state (which is a superposition of the two states with eigenvalue 
e?^) the x-component contributions will cancel, such that is given simply by the square of the z-component of the 
spin matrix expectation value on both sides of the QPT. This quantity is less than one throughout the broken phase, 
but approaches one at the QPT as can be seen immediately from Eqs. (24) and (25). As an illustration of this, we 
plot for the central-spin in Fig. |4]as a function oi X = fl ~ uJz- When A crosses 7a; = 1 in this diagram, we cross 
the QPT. From the numerics, we see that as the number of background spins is increased, the value approaches 
the ideal behavior in the thermodynamic limit. 



D. Transition between overlapping and non-overlapping energy surfaces 

For a certain range of parameter space the two energy surfaces may overlap along the energy axis, as we have 
already seen in Figs. [2](b, d). This also leads to a type of degeneracy in the system as different states on different 
sheets have equal energy. 

Note that e"^ forms the global maximum on the lower sheet and forms the global minimum on the upper sheet 
(this statement is true on both sides of the QPT). Hence at the point in parameter space where the transition to 
overlapping sheets occurs, these two points represent the first overlapping energy values. Therefore, making use of 
Eqs. ([l8| to set e*]. = immediately yields = as the transition point. Notice that these points must meet on 
the origin of the energy axis. 

Having obtained the condition = as the point where the transition to overlapping energy sheets occurs, we 
now divide our phase diagram for quadrant I into four distinct regions, as in Fig. [5j Here we have Regions 1 and 
2 as the broken phase region and 3 and 4 as the symmetric phase. The curve (w^ = acts as the dividing line 

between these regions. These two portions of the phase diagram are then further subdivided according to Regions 1 
and 3 as those without overlap and Regions 2 and 4 as those with overlapping energy sheets, with the line cj^ = 
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FIG. 5: (Color online) Here we show the phase diagram for Quadrant I from Fig. [T] in greater detail. We show the curve 

2 — 
tJz = ^ that denotes the quantum phase transition, as well as the line tj^ = which further divides the phase diagram into the 

odd-numbered regions (1, 3) with non-overlapping energy sheets and the even-numbered regions (2, 4) with overlapping sheets. 

We have chosen the value = 1. In the four surrounding panels we plot the density of states po against eo for representative 

values in each of the four regions. We choose the values as: Region 1 (f2 = 3/4, = 2); Region 2 {Q. = 2,liJz = 7/10); Region 

3 (SI = 1/3, = 1/2); Region 4 (f2 = 1/2, u)^ = 1/3), with 7a; = 1 in each case. 



acting as the dividing line in the latter case. 



IV. BEYOND MEAN FIELD THEORY: SPECTRA IN THE MAJORANA REPRESENTATION 



In this section we will relax the MFT approximation from the previous section and explicitly write the form of the 
Schrodinger equations in the Majorana representation as two de-coupled second-order differential equations. We will 
present numerical solutions to these equations on the Majorana sphere. We will then introduce the function G^'°'{a), 
as the logarithmic derivative of the Majorana polynomials <f>s,a{oi)- The roots of the Majorana polynomials will be 
transferred to the poles of the function, G"*'"(a), while the second-order equations will be transformed to first-order 
equations. We will then obtain an expansion for G^'°'{a) that allows us to obtain the integrated density of states in 
the thermodynamic limit in Sec. |Vj 



A. Schrodiner equation in terms of Majorana polynomials 



Having obtained the matrix form of the Hamiltonian ( 10 1 in the symmetrized basis, we can immediately re- write 
the Schrodinger equation in terms of the Majorana polynomials 



H(j){a) = e(j){a) 



(26) 
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in the form of two coupled first-order equations A^<f>s = {i^z/'^)<t'a and A-\-(j)a = {ujz/'i)4'si in which 

^ \ 2S ' 2S 

We can immediately de-couple these equations to obtain the second-order differential equations 

,2 



(27) 



4 



(28) 



in symbolic form. Now, making use of the expressions Q for the S operators in the coherent states representation 
and the definitions of A± given in Eq. ( 27 ) we can obtain the explicit form of these second order equations as 



5. + Po''"(a) 



(2S')2 ' 25- 
with the coefficients at each order given by 

((1 - 2S)a^ " 1) 7^ ± 2a07^ -\- 2S{2e - w + n){2e + uj + Q) 



(29) 



P 



ps,a _ 
^1 — 



8S 



{2S - l)a (a2 - 1) 7^ =F (a^ + l) ^Ix + 2an{n - 25(2e + n)) 



4S 



(30) 



Note that the symmetric and anti-symmetric expressions differ only with respect to terms that have no S dependence. 
Hence, in the thermodynamic limit these terms drop out and the equations for the two sectors will precisely agree. 
This implies that in the thermodynamic limit 0s and (j)a have the same functional dependence (hence also that tp^ and 
ip- must have the same form up to a constant in this limit). Therefore, we introduce the notation = lim5_>oo (t>^ — 
lims_>.oo (j)""- In Sec. (IVC) below we will explicitly consider this limit, after we have introduced the G function. 



B. Numerical plots of solutions on the Majorana sphere 

Here we plot some numerical results to illustrate our method and to gain insight into the physics of the system. 
This will also facilitate our discussion of the integrated density of states calculation in Sec. |V]and our analysis of the 
expectation value of system observables in Sec. |VI[ To conduct our simulations for a given set of parameter values, 
first we find the form of our (two) Majorana polynomials ([s]) by numerically solving the Mermin model Hamiltonian 
([T]). Then we project the roots of the polynomials onto the Majorana sphere using an inverse stereographic projection 
of the form 

-, / a + a \ 
Tria)^-—- {a-a)/i . (31) 
1 + V aa-l j 

One interesting feature of the Majorana sphere is that the arrangement of the roots projected onto the sphere reflects 
certain elements of the classical energy surface. In the present case, we will see parallels between the Majorana sphere 
projection and the classical energy surface diagrams presented in Fig. [2] To illustrate this point, in certain cases 
below we also project the level energy surfaces (the black curves in Fig. ^ onto the sphere. These are also visible as 
black curves in Figs. |6]-[9j 

In Fig. [6] we show projections of the numerically obtained roots for the low-excitation states near the ground state 
for values representative of Region 1. We use here the values lo^ = 2.5, = 0.5 and 7a; = 1 with the number of 
background spins N = 40. The roots of the 0s (a) polynomial are shown in blue and those of the 0a (o^) polynomial 
are shown in red. The first sphere on the far left represents the ground state, while the next three represent the first 
three excited states, respectively. 

In the case of the ground state at energy note that all of the poles line up along the imaginary axis in the upper 
hemisphere (in the complex plane, these poles extend to infinity along the positive and negative imaginary axis). As 
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FIG. 6: (Color online) Projection of the roots of the Majorana polynomial for low-excitation states in Region 1. The projections 
of the solutions to the 4>s{c() sector are shown in blue and those of the (j}a{a) sector are shown in red. Here we have chosen the 
numerical values tOz = 2.5, — 0.5 and = I (Region 1) with the number of background spins A'' — 40. From left to right, 
we have the k = 1 state (ground state), followed by A; = 2, 3, 4. The poles on the lower hemisphere line up along a straight-line 
contour on the real axis as we move from left to right. In each case, this contour is surrounded by the projected level energy 
surface, shown as the black circular curve. 



FIG. 7: (Color online) Projection of the roots of the Majorana polynomial for low-excitation states in Region 3. The projections 
of the solutions to the (j>s{a) sector are shown in blue and those of the (j}a{a) sector are shown in red. Here we have chosen the 
numerical values ujz = 0.95, SI — 0.5 and 7^, = 1 (region 3) with the number of background spins A*' = 40. From left to right, 
we have the fc = 1 state (ground state), followed by A; = 3, 5, 7 (we present odd numbers here to avoid a minor technicality). 
The two distinct level energy surfaces are shown in the A; = 3, 5, 7 spheres as the black curves. 



we move rightward we see that at each consecutive excited state a pair of poles vanishes from the upper hemisphere 
and appears at or in the vicinity of the origin (lower hemisphere). As they accumulate, we see that these poles are 
lining up on the real axis on the lower hemisphere as we move to higher excitation states. This process continues until 
all of the poles appear on the real axis, which occurs when we reach the upper edge of the lower energy sheet. 

This is useful for us because we can use the number of poles on a given (straight-line) curve as a way to index the 
polynomial solutions. Further on, we will define a complex integration contour around the poles on the real axis that 
allows us to 'count' which state we are in, and therefore we can calculate the integrated density of states as a function 
of energy e. As we consider other phase regions or different cases (such as overlapping sheets, etc.) we will find that 
this contour where the poles line up changes, however, some such contour(s) exists (exist) for all parameter values. It 
is this fact that allows us to obtain an analytic expression for the integrated density of states in each case. 

For the present case of Fig. [6j note further that the poles appearing in the lower hemisphere are enveloped by the 
level energy surface (black curve in spheres 2, 3 and 4). Here it is useful to think of moving rightwards among the 
spheres in Fig. |6]as moving upwards on the classical energy surface in Fig. [2|^a). Hence we may think of the poles on 
the lower hemisphere as 'nodes' that are expanding to fill the level surfaces, which themselves expand with increasing 
energy. 

In Fig. [7] we have given projections on the Majorana sphere for low-excitation states in Region 3 for the critical 
regime of the lower energy sheet. From left to right in this case we have the (broken phase) ground state at e5^, then 
the third, fifth and seventh excited states, respectively. In this case, the poles line up on two distinct contours along 
the real axis, but separated from the origin. Indeed these two contours are associated with the two potential wells on 
the lower energy sheet that we encountered previously in Fig. 2jc). These potential wells are represented in Fig. [t] 
by the two circular projected level energy surfaces enveloping the poles along the two straight-line contours. 

In Fig. [8]we present the Majorana sphere projections also for Region 3, but in this case we are passing the spectrum 
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FIG. 8: (Color online) Projection of the roots of the Majorana polynomial for excitation states near the saddle point in Region 
3. The projections of the solutions to the 0s (a) sector are shown in blue and those of the (j}a{c() sector are shown in red. Here 
we have chosen the numerical values uj^ = 0.95, Q = 0.5 and 7^ = 1 (region 3) with the number of background spins A'^ — 40. 
From left to right, we have the k — 9, 11, 17, 19 states. The level energy surface is shown in each sphere as the black curve. 




FIG. 9: (Color online) Projection of the roots of the Majorana polynomial for excitation states in the portion of the spectrum 
with overlapping sheets in Region 2. The projections of the solutions to the <^s(a) sector are shown in blue and those of the 
4>a{ci) sector are shown in red. Here we have chosen the numerical values ujz = 7/10, Q = 2 and yx ~ i (region 2) with the 
number of background spins A'^ = 40. From left to right, we have the k = 24, 25, 26, 27 states. Notice the alternating behavior 
of the "halo" shape, in the first sphere appearing on the lower hemisphere, then on the upper hemisphere in the second, and 
so forth. 

saddle-point at ( see Fig. ). Here it is clear that the two separate contours considered in Fig. [7] merge above 
the saddle-point to form a single contour. For the two spheres on the left of Fig. [8] we see that the two circular level 
surfaces have come together to form the united 'figure 8' shape from Fig. ^c). 

Finally in Fig. [9] we plot four spheres as an example of the behavior in the overlapping portion of the spectrum for 
Region 2 (here we have chosen to illustrate the k = 24, 25, 26, 27 excitation states). We have repressed the level energy 
surface curves in the two leftmost spheres. In this case, we see that as we move to higher excitation states, the "halo" 
shape appearing on the lower hemisphere in the leftmost sphere will appear alternately on the upper hemisphere 
(second sphere from the left) then return to the lower hemisphere (third sphere from the left) and so on. 

Again we can make sense of this behavior in terms of the level energy surfaces. Returning to the overlapping portion 
of the diagram in Fig. ^h) we see that the level energy surfaces appear as two circles: one small, inner circle (on the 
upper sheet) and one larger, outer circle (on the lower sheet). These represent two surfaces at the same energy that 
are not related by symmetry, similar to the previous studies of the LGM model (see the lower two spheres (e~, e^) in 
Fig. 4 of tl5J. 

These two surfaces are projected onto the two rightmost spheres in Fig. [9) Here we see that the halo shape seems 
to be alternating between the two sheets, lying more closely to the level surface in the lower hemisphere (upper sheet) 
in the third sphere from the left, while lying more closely to the level surface in the upper hemisphere (lower sheet) 
in the rightmost sphere. We will see that this alternating behavior again occurs in terms of the system observables in 
Sec. \VJ\ in particular see Figs. [T2| and [T4| . 
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C. G function and expansion in terms of particle number 

Let us now introduce the logarithmic derivative functions of our symmetric and anti-symmetric sector wave functions 



as 



25* 



fe=i 



„(s,a) ■ 



(32) 



where a^jf'"''' are respectively the roots of the symmetric and anti-symmetric polynomials. According to this transfor- 
mation, the roots of the Majorana polynomials have been transferred to poles of the functions G"''''(a). While G*'°(q!) 
is not exactly a Green's function, its roots serve to keep track of which state the system is in; hence, we can use it to 
calculate the integrated density of states. 



Introducing the G"*'"(q;) functions allows us to re-write the second order differential equations (29) as the first order 
equations 



Then, in order to take the thermodynamic limit, we can expand G^'" and e in the form of 



Si 



E 



(33) 



(34) 



With these expansions equation (32 1 reduces to a quadratic equation for Gq at the thermodynamic limit S* — >■ oo. 
These equations take the form 



Pr(a)(Go(a))' + Pr(a)Go(a) + Pg^ia) = 0. 
with the coefficients here obtained from Eqs. (pOl) as 



Po°^(a) = \{4e' + 4ne-u' + n'-a'^l) 
^r(«) = ^« 7^-2^^(26 + 17)) 

P^ia) = P^H = -f («' - PlKa^ Pi) 



where 



/3± 



(35) 

(36) 
(37) 
(38) 

(39) 



We emphasize that while we are still working in the thermodynamic limit here, we have relaxed the MFT assumption 
from Sec. |III[ hence we are working at a first-order level in the immediate approximation scheme. 



D. Solution for Go at the thermodynamic limit 

We can explicitly write the solution of the quadratic equation ( [35| for Go as 

a [217(260 + ^2)^2] ± 2 VO(a; eo) 

^0 = AP^) 

in which Q{a\ e) is itself a quadratic expression in given by 

with 

r_ (6) = (26 - - 17) (2e -f - 17)7^ 



(40) 

(41) 
(42) 
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The quadratic roots of Q{a; e) are given by 



7.(62- 6-2) ±a.,f]^(6^^)2 

7.(e-e5?)(e-e^) 



(43) 



in which 



472 ^ ^ 

(we repress the e dependence in r±(e) to avoid too heavy of notation). Recall that and are given in Eqs. im 
and (22), respectively. Since Q(a;eo) is under the radical in Eq. (40l, these roots describe two e-dependent b ranch 
cuts in the complex a plane; [— r_|_,r_|_] and [— r_,r_]. These cuts are precisely the contours discussed in Sec. IV B 
along which the roots will organize themselves with increasing e |21) . 

V. ANALYTIC CALCULATIONS FOR THE DENSITY OF STATES 

We have now developed the theoretical formalism that will allow us to perform the calculation of the integrated 
density of states, where the density of states itself is given by 



p(e) = ^5(6-e„) (45) 



for the Mermin model in the thermodynamic limit, where e„ are the system energy eigenvalues. This calculation is 
the primary goal of the present work. 

In order to facilitate our calculation, we define the density of poles 1(e) along the contour Ci that follows one of 
the branch cuts in the complex plane as 



I->^{e)^^^^[ G'^^ia-tjda (46) 
hi 



2S 2iTr 



in which p is the total number of poles along Ci and Ci is a contour that surrounds Ci. Ci itself follows one of the 
branch cuts in the complex plane, although this contour will change for different portions of the spectrum. 

Essentially what we have done is that we have transferred the zeroes of the Majorana polynomials into the poles 
of the functions G"''''(a; e). For given parameter values fi, uj^ and 7., as we increase e the positions of the poles will 
change as we saw in Sec. |IVB[ For example, as we increase the energy from the ground state in Region 1, these poles 
will gradually accumulate on the real axis in the lower hemisphere as in Fig. [6] Hence, by choosing our contour Ci 
along the the branch cut on the real axis from — to r_|_, the quantity in Eq. ( |46[ ) measures the accumulation of 
poles with increasing e. Since the pole density I^'°'{e) essentially counts the poles as we move from one state to the 
next, this quantity can be related to the integrated density of states 

AA(e) ^ r p[e')de' (47) 
J —00 

as in the following section. 

A. Explicit integral form of integrated density of states 

In the thermodynamic limit 5* — )■ 00 it can be shown that the integrated density of states A/'(e) can be written as 

m 

Afoieo) = lim Af'--{e) = lim ^'"(e) = lo(eo) (48) 

5— >-oo 5—^00 



in which eg is the first-order term for the energy defined in Eq. (34) and 



^o(eo) = ^/ da[G+{a;eo)-Goia;eo)]. (49) 
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The analytic portion of the integration wih cancel out, hence we have 

This is our standard integral form for the integrated density of states. In each case below, we can transform this 
integral in terms of the complete elliptic integral of the first kind 

K{k) = C '^'^ = (51) 

Jo ^{l-x-^){l-kx-^) 

and the complete elliptic integral of the third kind 

n(7i,fc)= / ^. (52) 

7o (l-na;2)y(l-a;2)(l-fca;2) 



B. Analytic form for integrated density of states by region 



Below we present the results of the integration given in Eq. (50) for each phase region in Fig. [5) We find that 
there are five different zones (a-e) corresponding to five different types of behaviors in the spectrum. Not all of these 
zone types actually occur in each phase region; only Region 4 is comprised of all five. These zones are shown in terms 
of the density of states for representative cases for all four regions in Fig. ([5]). In the Appendix |A] we show how to 
obtain the results presented in zone (d) by an analytic transformation of the integrand given in Eq. ( 50 1 , which is 
representative of the general method. 



1. Region 1 (symmetric, non- overlapping sheets) 

Let us first consider the case of Region 1 with non-degenerate ground state and non-overlapping energy surfaces. 
This is given by lo^ > 1x1^ ^'^'^ > $7 in Fig. [5j We find that the spectrum is comprised of zone (b) behavior on 
the lower sheet and zone (d) behavior on the upper sheet. 

(b) — '^'+^ < eo < (Lower sheet). The integrated density of states in this zone can be written in one of two 

ways, for example as 



<'^eo) = A^o(eo) 



(53) 



with 



No{e) 



aiR 



!± 

PI 



(54) 



(we remind the reader that there is an implicit dependence on e in r±(e) here). We have also defined 



ai,2 



(55) 



(We could also write the result as A/'g' (eo) = 5 " No{—eo) in which we have related the answer to that in zone (d) 
by symmetry). 



(d) ^ 



< eo < 



(Upper sheet). For the upper sheet, we can draw our contour Ci along the imaginary axis. 



The integrated DOS is given by 



(56) 



We show in App. [A|how to obtain this result by analytic manipulation of the integral in Eq. ( [50| . 

We plot our analytic results for Ao^(eo) for the full spectrum as the orange curve in Fig. 10 'a) against our numerical 
results for A/'^(e) (blue dots). We show the results for Regions 2, 3 and 4 in the panels labeled (b), (c) and (d), 
respectively. 
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[n=0.75; 6^=2.0; 7;,=!. 0] 



[n=2.0; 6^=0.70; 7;,=!. 0] 




(a) 



1.0 
0.8 



-1.0 -0.5 0.0 0.5 1.0 
[n=0.33; 6^=0.50; 7;,= 1.0] 




-0.6 -0.4-0.2 0.0 0.2 0.4 0.6 
(C) So 




(b) 



1.0 -0.5 0.0 0.5 1.0 

So 



[ n=0 .50; 6^=0.33; 7;,= 1.0] 




-0.6 -0.4-0.2 0.0 0.2 0.4 0.6 

(d) So 



FIG. 10: (Color online) Analytic (orange curve) and numerical (blue dots) integrated density of states in each region, with the 
parameter values as given above each plot (these are the same values as those used in Fig. [5] in each case), (a) Region 1, (b) 
Region 2, (c) Region 3, and (d) Region 4. 



2. Region 2 (symmetric, overlapping sheets) 

Now we consider the case of Region 2 with non-degenerate ground state and overlapping energy surfaces. This is 
given by lu^ > Ix/^ < in Fig. [s] The results are similar to that for Region 1, but here there occurs zone 

(c) behavior in the overlapping region of the spectrum. 

(b) — "'^^ < eo < (Non-overlapping portion of lower sheet). For this portion of the spectrum in Region 2, 

we have simply A/'o''^(eo) = ^^'^{eQ). 



(c) 



J, -a 



< eo < 



(Overlapping portion). For the overlapping portion of the spectrum we have 



(57) 



where we give N'^''^{e) immediately below, 
(d) 



< eo < "'2 (Non-overlapping portion of upper sheet). In the case of the upper sheet we have the 



results for A/'Q'^(eo) as in zone 1(d). 
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3. Region 3 (broken phase, non-overlapping sheets) 



Now we consider the case of region 3 with twice degenerate ground state (broken phase) and non-overlapping energy 
surfaces. This is given by utz < 7^/^^ and Wz > $7 in Fig. [sj Here the lower sheet is broken into zones (a) and (b), 
while the upper sheet is broken into (d) and (e). 



(a) 



■ll) 



< eo < 



(Critical portion of lower sheet). Throughout the degenerate portion 



of the lower sheet the contour Ci lies along two portions of the real axis as shown in Fig. [7] This results from the two 
branch cuts overlapping. In this case r_|_ > r_, so that the contour Ci extends from r_|_ to r_ then from — r_ to 



+ ■ 



Although in practice, it is easier to work out the integration in terms of the second contour that extends from infinity 
back to the origin along the positive imaginary axis before returning to infinity along the negative imaginary axis. 
We can obtain the integrated density of states in the critical region as 



A/'o'^(eo) 



+ iVo(-eo) - A^i(-eo) 



in which Afo(£o) is given in Eq. (54) and 



7Vi(e) 



_zyrl(^ 



- rir^ K 



(58) 



(59) 



We have obtained this result from symmetry in terms of that presented in zone 3(e) below, which can itself be obtained 
by following the method outlined in Appendix |A 3| 

At the precise value eo = — (w^ + ri)/2 that gives the boundary between the critical and non-critical parts of the 
spectrum on the lower sheet (i.e. the DOS singularity), we have r_ = and therefore the integrated DOS given in 
the form of Eq. (A4) reduces to elementary integrals. When computed, we find 



3 bJz + ^, 



I3l-r\ tan"^ 



^fil-r\ tan-i 




(60) 



(b) — < eo < (Non-critical portion of lower sheet). For the non-critical portion of the lower sheet in 

Region 3, we have the usual result for A/'o'''(eo) the same as in Eq. (53). 



(d) 



< eo < 



(Non-critical portion of upper sheet). In the case of the upper sheet wc have the typical 



behavior for M^'^ito) given in Eq. (56). 

At the boundary cq = (w^ -I- r2)/2 between the critical and non-critical parts of the spectrum on the upper sheet, 
the integrated DOS expression again simplifies. In this case wc have the result 



AAo^(^)-l-AA3(-^ 



(61) 



in terms of Eq. (60). 

(e) < eg < — 2^-\/(a;2 -^^'^(^'^ +7^) (Critical portion of upper sheet). In this case, we have the result 



A/'o''(eo)= 2-^o(eo) + A^i(eo) 



(62) 



as outlined in the Appendix |A 3 



J^. Region 4 (broken phase, overlapping sheets) 



Now we consider the case of Region 4 with degenerate ground state and overlapping energy surfaces. This is given 
hy ujz < Jx/^ s-iid a;^ < 51 in Fig. [5] Here the spectrum consists of zones (a) and (b) on the lower sheet, then an 
overlapping region (c) and then the remaining portion of the upper sheet consists of zones (d) and (e). 

(a) —i^^(uj1 -|-72)(r22 4-7^) < eo < — "'^^^ (Critical portion of lower sheet). Here the spectrum imitates exactly 



that from zone 3(a) with '^(eo) — A/'q 



At the critical boundary cq 



■^'^(eo) as given in Eq. (58) 
f2)/2 we again have N^{~ 



2 ' 



2 



) as given in Eq. ( 60 ) . 
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FIG. 11: (Color online) Comparison of expectation values for the observables appearing in the Hamiltonian Eq. ([T| for Region 1. 
Respectively by column we have plotted {(Jz), (Sz) and (a^Sx) across the range of both energy sheets for the system parameter 
values ljJz = 2, Q — 3/4 and yx = 1- The green curves represent the results from the Hellmann-Feynman theorem and the blue 
dots are numerical results computed assuming A'' = 40 as the number of background spins. The red dashed lines indicate the 
position of the energy ^b± as given in Eq. (66 1. 



(b) 
T TV 
(c) 



(Non-overlapping portion of non-critical lower sheet). Here again we have the usual result 



< eo < 2 

E 

"^"^^ (Overlapping portion). For the overlapping spectrum we have 



for A/'o''^(eo) as given in Eq. (53) 



< £0 < 



A/'o'^M=2-(A/'o''M+<-'^(eo)) 



(63) 



similar to the result given in zone 2(c) in Eq. (57). 

(d) < eo < (Non -overlapping portion of non-critical upper sheet). In the case of the non-overlapping, 
non-critical part of the upper sheet we have again the typical behavior for M^'^ito) given in Eq. (56). 

At the upper critical boundary eo = (w^ -f ri)/2 we again have ( "^^^ ) = 1 - Mq{-'^^^^) similar to Eq. (60) in 
Region 3. 

(e) < eo < — 2^-\/(w2 +72)(ri2 +^^2^ (Critical portion of upper sheet). Here we again have M^''^((.q) = 
A/'o'''(eo) as given in Eq. (62). 



VI. EXPECTATION VALUES OF SYSTEM OBSERVABLES 



In this section we present the expectation values for several key observables in the Mermin model. For the most 
important observables (those appearing in the Hamiltonian Eq. ([T])) we can easily obtain expectation values from our 
analytic results by applying the Hellmann-Feynman theorem. The Hellmann-Feynman theorem relates the expectation 
values of observables to partial derivatives with respect to the system parameters as they appear in the energy 
eigenvalues. In our case we easily obtain, for example, 

, > „ 9e „ 9e dM 

where we have introduced partial derivatives in terms of the integrated density of states M since our expressions for 
this quantity given in Sec. |VB| are too complicated to invert analytically for e. 



In Fig. 11 we plot the expectation values for (ctz), {S^) and {(JxSx), respectively across the full spectrum for a 
representative case in Region 1 (symmetric case with non-overlapping sheets). We have also included the numerical 
results, calculated directly from the Majorana polynomials using a system with = 40 background spins. In the 
lower left panel we have plotted the expectation value for the central spin {a^) across the energy range of the lower 
sheet. We see that at energy values e near the ground state, this spin points upward, consistent with the picture 
presented in Fig. |3][a). Meanwhile, in the middle lower diagram the background spin {Sz) also points up at energies 
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<a> <S^>/S <a^S^>/S 




••: >j -i.o i •• >j 1 -vvvv^- , 

-1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 -1.0 -0.5 0.0 0.5 1.0 



FIG. 12: (Color online) Comparison of expectation values for the observables appearing in the Hamiltonian Eq. ([T]) for Region 2. 
Respectively by column we have plotted (cr^), {S^} and (a^S^) across the spectrum for the system parameter values lj^ = 7/10, 
n = 2 and = 1- The green curves represent the results from the Hellmann-Feynman theorem and the blue dots are numerical 
results computed assuming = 40 as the number of background spins. 



near the ground state. As we consider higher excitation states, we see that the z component of both spins decreases as 
the effect of the coupling tends to pull the two spins towards anti-alignment on the s-axis. However, the central spin 
reaches a minimum value partway across the range of the lower sheet and begins to re-align with the positive z-axis 
(the approximate energy value es- where this occurs is labeled with a dashed red line). Meanwhile, the z component 
for the background spin continues leveling out, until it actually becomes negative in the upper half of the lower sheet. 
Finally, at the top of the lower sheet the two spins are anti-aligned in the z-axis. This is consistent with the picture 
presented in the lower right-hand panel for {uxSx)- Here we see that the anti- alignment of the x-component of the 
spins is zero at the ground state, reaches a maximum partway across the energy sheet at es- then returns to zero at 
the top of the lower energy sheet. 

A similar behavior occurs in the range of the upper sheet, only with the central spin having a negative instead of a 
positive alignment on the z axis. Also {ct^Sx) < in the upper sheet, which indicates that in this case alignment on 
the a;-axis occurs between the two spins, instead of anti-alignment. 

In theory, we could determine the spectruml extremal points that give the values of e±B using our expressions for 



the integrated density of states given in Eqs. (531 and (56), but the expressions are too complicated to invert them to 



explicitly obtain the spectrum. However, we can still determine the approximate analytic position of €±b using our 



mean field theory expression for the energy given in Eq. (11). We determine the extremal points in the expectation 
value for according to 

This condition yields the 'breaking point' values for a, a according to as — olb = il- Plugging these values back 



into Eq. (11) yields the breaking point energies as 



eB.^±^^^. (66) 

Note that there is no dependence on il in this expression, only the coupling to the environment affects the position 
of these extremal points at this order of approximation. These points are labeled by the red dashed lines appearing 
in Figs. [TT| 

In Fig. |12|we plot these same three expectation values for a case representative of Region 2, including numerical 
results. Here each plot can be broken up into three sections along the lines of Sec. |VB 1| We see immediately that 
the numerical results give us two distinct "bands," even in the overlapping region at the horizontal center of each 
plot. This can be viewed in relation to Fig. [9| We have seen in this figure that increasing energies give an alternating 
'halo' pattern of zeros on the Majorana sphere in the overlapping region; similarly, we see in Fig. [12] that increasing 
energies give states with expectation values that alternate between the two sheets in this region. 

However, the analytic results obtained from application of the Hellmann-Feynman theorem (green curve) give an 
averaged value in the overlapping region, hence we would lose some information if we were to rely solely on Hellmann- 
Feynman. This behavior is analogous to that occuring in Region IV of the LMG model as reported in 14 . 



We plot the same set of expectation values for a representative case of Region 3 in Fig. 13 for oj^ = 1/2, VL = 1/3 and 
7a: = 1. We see that these expectation values are sharply peaked at the critical point at eL,e^ ~ T5/12 ~ =fO. 416667 
between the broken and symmetric portions of the spectrum, as we should expect. 

Finally we plot these expectation values for the system observables in Region 4 in Fig. |14[ Not surprisingly, we see 
a combination of the behaviors previously occurring in Regions 2 and 3. 
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FIG. 13: (Color online) Comparison of expectation values for the observables appearing in the Hamiltonian Eq. (|T]) for Region 
3. Respectively by column we have plotted {(Jz), {Sz) and (a^Sx) across the range of energies for both sheets, as labeled, for the 
system parameter values uuz = 1/2, il = 1/3 and 7a; = 1. The green curves represent the results from the Hellmann-Feynman 
theorem and the blue dots are numerical results computed assuming A*' = 40 as the number of background spins. 
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FIG. 14: (Color online) Comparison of expectation values for the observables appearing in the Hamiltonian Eq. ([T]) for Region 4. 
Respectively by column we have plotted (ctz), (Sz) and (a^Sx) across the spectrum for the system parameter values Uz = 1/3, 
Q = 1/2 and 73; = 1. The green curves represent the results from the Hellmann-Feynman theorem and the blue dots are 
numerical results computed assuming A*' = 40 as the number of background spins. 



VII. CONCLUSION 



In this study we have analyzed Mermin's special case of the central-spin model with symmetric coupling to environ- 
mental modes and monochromatic bath in the framework of a spin-coherent states formalism. By applying a mean 
field theory approximation, we were able to obtain the classical energy surfaces for the double spectrum as presented 
in Eq. (11) and plotted for representative cases in Fig. [2j our results here are in agreement with studies by other 
researchers [SHZ]. We also obtained a complete phase picture for the model as presented in Figs, [l] [sj including the 
quantum phase transition jw^j = 7;?/|f^| between Regions 1, 2 (symmetric phase) and Regions 3, 4 (broken phase) 
as well as the crossover point \ujz\ ~ \^\ for overlapping energy sheets between Regions 1, 3 (non-overlapping) and 
Regions 2, 4 (overlapping). 

Moving beyond the mean field theory approach, from the time-independent Schrodinger equation we were able to 
obtain the Gq function in the thermodynamic limit in Eq. ( 40 ) . Relying on the integration contours we discovered 
from Figs. [6] -[9) we were able to write a closed form expression for the integrated density of states (for each phase 
region) in terms of an integration of the non-analytic part of the Gq function over one of these contours. This integral 
is presented in generic form in Eq. (50 1. We used this result in Sec. VB to work out the integrated density of states 

Some details of how to obtain these 
Finally in Sec 



VI 



in analytic form for every zone of the spectrum in each of the four phase regions 
expressions by analytic manipulations of the integral in Eq. ( 50 ) are presented in the Appendix 
we plotted numerical results for the observables appearing in the Hamiltonian for representative cases in each of the 
four regions and applied the Hellmann-Feynman theorem to obtain an accurate comparison from our analytic results 
for the integrated density of states. 
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We envision several avenues for future study on the Mermin model as well as related models. For example, in the 
Mermin model it might be interesting to consider in greater detail the 1/S corrections to the G^^" function in Eq. 
(34), as it is at first order in l/S" that the symmetric and anti-symmetric functions can be differentiated from one 
another according to Eqs. (30). While these calculations would be involved, they would allow us to determine the 
spectruml differences between the two sectors. 

Dynamical studies of the Mermin model may also provide the opportunity to study the evolution of a two-level 
system in the presence of a non-Markovian bath. In particular, a semi-classical analysis building from the formalism 
presented in this paper may prove fruitful. Other approaches to the dynamics of the Mermin model have already been 
explored in Refs. [5] and The presence of non-Markovian effects has been demonstrated in [S]. 

One further well-posed problem would be to study the generalization of the present model to a central system with 
s = 1, or even a two qubit central system for which both the internal entanglement as well as that between the central 
system and the bath could be computed. 
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Appendix A: Explicit calculations related to the integrated density of states 
1. Transformation of integrated density of states in non-critical zone (d) 



In this section of the appendix we demonstrate how to manipulate the integral in Eq. ( 50 ) to obtain the result given 



in Eq. (56) for the (d) zone of the DOS. Setting aside an eo-dependent coefficient, the integral we should consider is 
given by 



;o (a2-/32)(a2_/3|) 
It is known that any integral over a rational function of a and 



da. 



y{a) = \J {a^ — r\){a'^ — r^) 



(Al) 



(A2) 



can be written as a combination of the three types of elliptic integrals in standard form |19j (in general {y{a))'^ should 
be any cubic or quartic in a for this statement to apply). 

To obtain the integral above in standard form, we assume that the integrand can be written as a sum of integrands 
in the form 



B 



C 



PI) y{a) (o^ - /3|)y(a) 



+ Dy{a) 



(A3) 



in which A, B, C and D are undetermined constants. After slight re-arrangement these terms would correlate to the 
integrands for the elliptic integral of the first kind, two possible elliptic integrals of the third kind and the elliptic 
integral of the second kind, respectively. To determine these unknown constants we multiply Eq. (A3) through by 
a factor of y{a){a'^ — l3i)(a'^ — /3|), yielding a quartic polynomial in from which we can determine the vmknowns 
A,B,C and D. We immediately obtain D — (no elliptic integral of the second kind) as well as A = 1, i? = /3^ai 



ia2, yielding the coefficients given in Eq. (55) 



and C = 

From this point, re-arranging each integrand slightly and performing a transformation of the integration variable 
according to q = oijr^ allows us to obtain the result in Eq. (56) in terms of elliptic integrals in standard form. 
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2. Integration form at the saddle-points in Regions 3 and 4 



In general, the integral form for the integrated DOS in Eq. (50) can be re- written after a simple fraction decompo- 
sition as 



1 



[A(A)-A(/32)] 



(A4) 



with the integral A(/3) given by 



m 



Ci 



a2 -^2 



da. 



(A5) 



In the special case eo = — (wz -I- fi)/2 in Regions 3 and 4 we have r_ = and the integral A(/3) reduces so that it can 
be computed through elementary techniques, yielding Eq. (60 1. 



3. Transformation of DOS integral in the critical zone (e) 



In this region, it is easiest to re- write the integral in terms of a contour that extends from infinity back to the origin 
along the positive imaginary axis and then back out to infinity on the negative imaginary axis. We can re-write the 
integration variable in Eq. (50) along the imaginary axis according to iq = a. Noting the symmetry properties of the 



resulting integral, it can be written as proportional to 



(g2 



I X 



-dq. 



(A6) 



We can break Yo{e) into two integrals, one over (0, r+) and a second over (r+, oo). The first integral (re-including all 
appropriate coefficients) yields Nq from Eq. (54). 

For the remaining integral on the interval (r+, oo) we perform a simple inversion of the integration variable according 
to — q. After some slight re-arranging of the resulting integrand, we obtain a new integral along (0, 1/?'+) that 
can again be calculated by following a method similar to that outlined in App. |A 1| above. After re-introducing 
coefficients, we obtain the result Ni(e) as given in Eq. (59) in the main text. 
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